%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%This m-file generates Tables 10, 11 and Figure 3
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all
clc
load BayesBeta_normal_3.mat

Age=Demo(:,1)';

K=3;
D=3;

ncomp3=Bayesbeta(:,1:K)';                        %K by S matrix
ind3=Bayesbeta(:,K+1:K+N)';                      %N by S matrix
pvec3=Bayesbeta(:,K+N+1:K+N+K)';                 %K by S matrix
pvec3mean=mean(pvec3(:,5000:S),2);
Beta_test=Bayesbeta(:,K+N+K+1:K+N+K+L*N);   
Beta3=reshape(Beta_test',L,N,S);                 %L by N by S matrix  
clear Beta_test
Beta3mean=reshape(mean(reshape(Beta3(:,:,5000:S),L*N,[]),2),L,N);
beta_test=Bayesbeta(:,K+N+K+L*N+1:K+N+K+L*N+L*D);
beta3=reshape(beta_test',D,L,S);                      %D by L by S matrix
clear beta_test
beta3mean=reshape(mean(reshape(beta3(:,:,5000:S),D*L,[]),2),D,L);
beta3std=reshape(std(reshape(beta3(:,:,5000:S),D*L,[]),0,2),D,L);
beta3p5=reshape(prctile(reshape(beta3(:,:,5000:S),D*L,[]),5,2),D,L);
beta3p95=reshape(prctile(reshape(beta3(:,:,5000:S),D*L,[]),95,2),D,L);
b_test=Bayesbeta(:,K+N+K+L*N+L*D+1:K+N+K+L*N+L*D+K*L);
b3=reshape(b_test',L,K,S);                             %L by K by S matrix
clear b_test
b3mean=reshape(mean(reshape(b3(:,:,5000:S),L*K,[]),2),L,K);
W_test=Bayesbeta(:,K+N+K+L*N+L*D+K*L+1:K+N+K+L*N+L*D+K*L+L*L*K);
W3=reshape(W_test',L,L,K,S);                          %L by L by K by S matrix
clear W_test
W3mean=reshape(mean(reshape(W3(:,:,:,5000:S),L*L*K,[]),2),L,L,K);

Beta3=Beta3(:,:,5001:S);
clear Bayesbeta



% The age distribution in HK (from Census 2006)
%A=9;                        %number of age group in Census (from age 15-59)
age_group(1)=223369+215868; %population size in age_group 1
age_group(2)=224834+244934;
age_group(3)=223704+278583;
age_group(4)=239021+310818;
age_group(5)=247397+329963;
age_group(6)=305447+366048;
age_group(7)=323497+336158;
%age_group(8)=264753+269380;
%age_group(9)=214652+207650;
A=size(age_group,2);


G=5;                        %number of age in an age group
for a=1:A
    for g=1:G
        age_prob((a-1)*G+g)=age_group(a)/(sum(age_group)*G);
    end
end

% The family income distribution in HK (frmo Census 2006)
FI=7;                       %number of family income group
income_group(1)=86736;
income_group(2)=118779+121605+146010+147081;
income_group(3)=339469+279217;
income_group(4)=225292+162783;
income_group(5)=221101;
income_group(6)=194723;
income_group(7)=183750;

for i=1:FI
    income_prob(i)=income_group(i)/sum(income_group);
end

% BT distribution in HK (from IP report 2008 of Dept of IP in HK)
BT_prob=.3;


% Do the simulations 100 times
SIM=100;

% Simulate 1000 individuals according to the demographics assumed
Sim=10000;

age=14+randsample(A*G,Sim,true,age_prob);
income=-1+randsample(FI,Sim,true,income_prob);
bt=-1+randsample(2,Sim,true,[1-BT_prob BT_prob]);

Demo_sim=[age bt income];

for ss=1:SIM



% Simulate the unobserved heterogeneity
ind_sim=randsample(K,Sim,true,pvec3mean');
for s=1:Sim
    test=10;
    beta_price_sim=Demo_sim(s,:)*beta3mean(:,4);
    while test+beta_price_sim>0
        rBeta_sim(s,:)=mvnrnd(b3mean(:,ind_sim(s))',W3mean(:,:,ind_sim(s)));
        test=rBeta_sim(s,4);
    end
end

temp_beta=mvnrnd(reshape(beta3mean,1,[]),diag((reshape(beta3std,1,[])).^2));
temp_beta=reshape(temp_beta,3,[]);
temp=Demo_sim*temp_beta+rBeta_sim;
%Beta_sim=Demo_sim*beta3mean;
Beta_sim(:,:,ss)=temp';

clear temp temp_beta
end
% Calculate the marginal cost using the original coefficient
elas3=elasticity(Beta3(:,Age<22,:),Scenario);
mc=Scenario(1,1)*(1+1/elas3(1,1));

Scenario_adult=Scenario;
Scenario_adult(1,1)=10;

% Calculate the change in shares and profit under different counterfactuals
test=counterfactual(Beta3, Scenario,mc, [1 1 1]);
share_all_available=test(:,1:4);
profit_all_available=test(:,5:8);
clear test
test=counterfactual(Beta3, Scenario,mc, [1 0 1]);
share_no_street=test(:,1:4);
profit_no_street=test(:,5:8);
clear test
test=counterfactual(Beta3, Scenario,mc, [1 1 0]);
share_no_internet=test(:,1:4);
profit_no_internet=test(:,5:8);
clear test
test=counterfactual(Beta3, Scenario,mc, [1 0 0]);
share_no_piracy=test(:,1:4);
profit_no_piracy=test(:,5:8);
clear test

test=counterfactual(Beta_sim(:,age'<22,:), Scenario,mc, [1 1 1]);
share_stu_all_available=test(:,1:4);
profit_stu_all_available=test(:,5:8);
clear test
test=counterfactual(Beta_sim(:,age'<22,:), Scenario,mc, [1 0 1]);
share_stu_no_street=test(:,1:4);
profit_stu_no_street=test(:,5:8);
clear test
test=counterfactual(Beta_sim(:,age'<22,:), Scenario,mc, [1 1 0]);
share_stu_no_internet=test(:,1:4);
profit_stu_no_internet=test(:,5:8);
clear test
test=counterfactual(Beta_sim(:,age'<22,:), Scenario,mc, [1 0 0]);
share_stu_no_piracy=test(:,1:4);
profit_stu_no_piracy=test(:,5:8);
clear test

test=counterfactual(Beta_sim(:,age'>=22,:), Scenario_adult,mc, [1 1 1]);
share_adult_all_available=test(:,1:4);
profit_adult_all_available=test(:,5:8);
clear test
test=counterfactual(Beta_sim(:,age'>=22,:), Scenario_adult,mc, [1 0 1]);
share_adult_no_street=test(:,1:4);
profit_adult_no_street=test(:,5:8);
clear test
test=counterfactual(Beta_sim(:,age'>=22,:), Scenario_adult,mc, [1 1 0]);
share_adult_no_internet=test(:,1:4);
profit_adult_no_internet=test(:,5:8);
clear test
test=counterfactual(Beta_sim(:,age'>=22,:), Scenario_adult,mc, [1 0 0]);
share_adult_no_piracy=test(:,1:4);
profit_adult_no_piracy=test(:,5:8);
clear test

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Changes in Estimated Demand and Profits When Street Piracy Not Available
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

demand_change_legal=(share_no_street(:,1)-share_all_available(:,1))./share_all_available(:,1);
demand_change_internet=(share_no_street(:,3)-share_all_available(:,3))./share_all_available(:,3);
demand_change_stu_legal=(share_stu_no_street(:,1)-share_stu_all_available(:,1))./share_stu_all_available(:,1);
demand_change_stu_internet=(share_stu_no_street(:,3)-share_stu_all_available(:,3))./share_stu_all_available(:,3);
demand_change_adult_legal=(share_adult_no_street(:,1)-share_adult_all_available(:,1))./share_adult_all_available(:,1);
demand_change_adult_internet=(share_adult_no_street(:,3)-share_adult_all_available(:,3))./share_adult_all_available(:,3);

profit_change_no_street=(profit_no_street(:,1)-profit_all_available(:,1));
profit_change_stu_no_street=(profit_stu_no_street(:,1)-profit_stu_all_available(:,1));
profit_change_adult_no_street=(profit_adult_no_street(:,1)-profit_adult_all_available(:,1));
profit_change_no_street_BSA=(Scenario(1,1)-mc)*share_all_available(:,2);
profit_change_stu_no_street_BSA=(Scenario(1,1)-mc)*share_stu_all_available(:,2);
profit_change_adult_no_street_BSA=(Scenario_adult(1,1)-mc)*share_adult_all_available(:,2);


Table_10_estimates = zeros(4,3);
Table_10_estimates(1,1)=mean(demand_change_legal);
Table_10_estimates(2,1)=mean(demand_change_internet);
Table_10_estimates(3,1)=mean(profit_change_no_street)*100;
Table_10_estimates(4,1)=mean(profit_change_no_street_BSA)*100;
Table_10_estimates(1,2)=mean(demand_change_stu_legal);
Table_10_estimates(2,2)=mean(demand_change_stu_internet);
Table_10_estimates(3,2)=mean(profit_change_stu_no_street)*100;
Table_10_estimates(4,2)=mean(profit_change_stu_no_street_BSA)*100;
Table_10_estimates(1,3)=mean(demand_change_adult_legal);
Table_10_estimates(2,3)=mean(demand_change_adult_internet);
Table_10_estimates(3,3)=mean(profit_change_adult_no_street)*100;
Table_10_estimates(4,3)=mean(profit_change_adult_no_street_BSA)*100

Table_10_prctile5 = zeros(4,3);
Table_10_prctile5(1,1)=prctile(demand_change_legal,5);
Table_10_prctile5(2,1)=prctile(demand_change_internet,5);
Table_10_prctile5(3,1)=prctile(profit_change_no_street,5)*100;
Table_10_prctile5(4,1)=prctile(profit_change_no_street_BSA,5)*100;
Table_10_prctile5(1,2)=prctile(demand_change_stu_legal,5);
Table_10_prctile5(2,2)=prctile(demand_change_stu_internet,5);
Table_10_prctile5(3,2)=prctile(profit_change_stu_no_street,5)*100;
Table_10_prctile5(4,2)=prctile(profit_change_stu_no_street_BSA,5)*100;
Table_10_prctile5(1,3)=prctile(demand_change_adult_legal,5);
Table_10_prctile5(2,3)=prctile(demand_change_adult_internet,5);
Table_10_prctile5(3,3)=prctile(profit_change_adult_no_street,5)*100;
Table_10_prctile5(4,3)=prctile(profit_change_adult_no_street_BSA,5)*100

Table_10_prctile95 = zeros(4,3);
Table_10_prctile95(1,1)=prctile(demand_change_legal,95);
Table_10_prctile95(2,1)=prctile(demand_change_internet,95);
Table_10_prctile95(3,1)=prctile(profit_change_no_street,95)*100;
Table_10_prctile95(4,1)=prctile(profit_change_no_street_BSA,95)*100;
Table_10_prctile95(1,2)=prctile(demand_change_stu_legal,95);
Table_10_prctile95(2,2)=prctile(demand_change_stu_internet,95);
Table_10_prctile95(3,2)=prctile(profit_change_stu_no_street,95)*100;
Table_10_prctile95(4,2)=prctile(profit_change_stu_no_street_BSA,95)*100;
Table_10_prctile95(1,3)=prctile(demand_change_adult_legal,95);
Table_10_prctile95(2,3)=prctile(demand_change_adult_internet,95);
Table_10_prctile95(3,3)=prctile(profit_change_adult_no_street,95)*100;
Table_10_prctile95(4,3)=prctile(profit_change_adult_no_street_BSA,95)*100



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Changes in Estimated Demand and Profits When All Piracy Not Available
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

profit_change_no_piracy=(profit_no_piracy(:,1)-profit_all_available(:,1));
profit_change_stu_no_piracy=(profit_stu_no_piracy(:,1)-profit_stu_all_available(:,1));
profit_change_adult_no_piracy=(profit_adult_no_piracy(:,1)-profit_adult_all_available(:,1));

profit_change_no_street_BSA=(Scenario(1,1)-mc)*share_all_available(:,2);
profit_change_no_internet_BSA=(Scenario(1,1)-mc)*share_all_available(:,3);
profit_change_no_piracy_BSA=profit_change_no_street_BSA+profit_change_no_internet_BSA;
profit_change_stu_no_street_BSA=(Scenario(1,1)-mc)*share_stu_all_available(:,2);
profit_change_stu_no_internet_BSA=(Scenario(1,1)-mc)*share_stu_all_available(:,3);
profit_change_stu_no_piracy_BSA=profit_change_stu_no_street_BSA+profit_change_stu_no_internet_BSA;
profit_change_adult_no_street_BSA=(Scenario_adult(1,1)-mc)*share_adult_all_available(:,2);
profit_change_adult_no_internet_BSA=(Scenario_adult(1,1)-mc)*share_adult_all_available(:,3);
profit_change_adult_no_piracy_BSA=profit_change_adult_no_street_BSA+profit_change_adult_no_internet_BSA;

Table_11_estimates = zeros(2,3);
Table_11_estimates(1,1) = mean(profit_change_no_piracy)*100;
Table_11_estimates(1,2) = mean(profit_change_stu_no_piracy)*100;
Table_11_estimates(1,3) = mean(profit_change_adult_no_piracy)*100;
Table_11_estimates(2,1) = mean(profit_change_no_piracy_BSA)*100;
Table_11_estimates(2,2) = mean(profit_change_stu_no_piracy_BSA)*100;
Table_11_estimates(2,3) = mean(profit_change_adult_no_piracy_BSA)*100

Table_11_prctile5 = zeros(2,3);
Table_11_prctile5(1,1) = prctile(profit_change_no_piracy,5)*100;
Table_11_prctile5(1,2) = prctile(profit_change_stu_no_piracy,5)*100;
Table_11_prctile5(1,3) = prctile(profit_change_adult_no_piracy,5)*100;
Table_11_prctile5(2,1) = prctile(profit_change_no_piracy_BSA,5)*100;
Table_11_prctile5(2,2) = prctile(profit_change_stu_no_piracy_BSA,5)*100;
Table_11_prctile5(2,3) = prctile(profit_change_adult_no_piracy_BSA,5)*100


Table_11_prctile95 = zeros(2,3);
Table_11_prctile95(1,1) = prctile(profit_change_no_piracy,95)*100;
Table_11_prctile95(1,2) = prctile(profit_change_stu_no_piracy,95)*100;
Table_11_prctile95(1,3) = prctile(profit_change_adult_no_piracy,95)*100;
Table_11_prctile95(2,1) = prctile(profit_change_no_piracy_BSA,95)*100;
Table_11_prctile95(2,2) = prctile(profit_change_stu_no_piracy_BSA,95)*100;
Table_11_prctile95(2,3) = prctile(profit_change_adult_no_piracy_BSA,95)*100


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Simulated Profit Under Different Piracy Prices (Figure 3)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
demand=zeros(6,6);
profit=zeros(6,6);
for i=0:1:5
    for j=0:1:5
        Scenario(2,1)=i;
        Scenario(3,2)=j;
        test=counterfactual(Beta3, Scenario,mc, [1 1 1]);
        demand(i+1,j+1)=mean(test(:,1));
        profit(i+1,j+1)=mean(test(:,5));
    end
end

profit=profit*100;


figure(1)
colormap gray
mesh(0:1:5,0:100:500,profit)
xlabel('Days of Downloading', 'fontsize',12)
ylabel('Price for Street-Pirated Office','fontsize',12)
zlabel('Expected Profit Per Student ($HK)','fontsize',12)
break